Draft version February 3, 2008 

Preprint typeset using I^T^X style cmulatcapj v. 6/22/04 



AN EVOLUTIONARY MODEL FOR SUBMILLIMETER GALAXIES 
Sukanya Chakrabarti 1,2 , Yeshe Fenner 1 , Lars Hernquist 1 , T.J. Cox 1 , Philip F. Hopkins 1 

Draft version February 3, 2008 

ABSTRACT 

We calculate multi-wavelength spectral energy distributions (SEDs) (spanning optical to millimeter 
wavelengths) from simulations of major galaxy mergers with black hole feedback which produce sub- 
millimeter bright galaxies (SMGs), using a self-consistent three-dimensional radiative transfer code. 
These calculations allow us to predict multiwavelength correlations for this important class of galaxies. 
We review star formation rates, the time evolution of the 850 fim fluxes, along with the time evolution 
of the Mbh — M s t al relation of the SMGs formed in the mergers. We reproduce correlations for local 
AGN observed in Spitzer Space Telescope's IRAC bands, and make definitive predictions for infrared 
X-ray correlations that should be testable by combining observations by Spitzer and the upcoming 
Herschel mission with X-ray surveys. To aid observational studies, we quantify the far-infrared X-ray 
correlations. Our dynamical approach allows us to directly correlate observed clustering in the data 
as seen in IRAC color-color plots with the relative amount of time the system spends in a region of 
color-color space. We also find that this clustering is positively correlated with the stars dominating in 
their contribution to the total bolometric luminosity. The merger simulations also allow us to directly 
correlate the 850 /im flux with the ratio of the black hole luminosity to the total luminosity, which 
is an inherent and testable feature of our model. We present photo albums spanning the lifetime of 
SMGs, from their infancy in the pre-merger phase to the final stage as an elliptical galaxy, as seen 
in the observed 3.6 /xm and 450 /im band to visually illustrate some of the morphological differences 
between mergers of differing orbital inclination and progenitor redshift. We compare our SEDs from 
the simulations to observations of SMGs and find good agreement. We find that SMGs are a broader 
class of systems than starbursts or quasars. We introduce a simple, heuristic classification scheme 
on the basis of the Lir/L x ratios of these galaxies, which may be interpreted qualitatively as an 
evolutionary scheme, as these galaxies evolve in Lir/L x while transiting from the pre-merger stage, 
through the quasar phase, to a merger remnant. 

Subject headings: galaxies: formation — galaxies: AGN — infrared: galaxies — radiative transfer — stars: 
formation 



1. INTRODUCTION 

The luminous (Lir, > 1 x 10 12 L Q ), high-redshift, 
dusty galaxies which were discovered in large num- 
bers by the Sub-millimeter Common-User Bolometer Ar- 
ray (SCUBA), now designated as sub-millimeter galaxies 
(SMGs), generate a significant fraction of the cosmic en- 
ergy output (Smail et al. 1997; Ivison et al. 1998; Blain 
et al. 2002). As such, SMGs are key cosmological play- 
ers, believed to be responsible for more than half of the 
star formation at z ~ 2 (Blain et al. 2002). 

These systems have been studied now in a diverse range 
of wavelengths, spanning the range from X-rays to radio 
wavelengths. UV/optical spectroscopic redshifts (Chap- 
man et al. 2003b, 2005, Swinbank et al. 2004) yield a 
median redshift of ~ 2 for this population. The presence 
of AGN in these systems has been confirmed in X-ray 
surveys (Alexander et al. 2005a, b), however the contri- 
bution of the AGN to the bolometric luminosity remains 
unclear. Borys et al.'s (2005) analysis of rest-frame X-ray 
and near-IR data suggests that SMGs fall nearly two or- 
ders of magnitude below the local M s t a r — -^bh relation, 
when Eddington-limited accretion is assumed to calcu- 
late the black hole masses. HST imaging indicates that 
these are large irregular galaxies (Chapman et al. 2003a), 
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with some preliminary evidence for extended starbursts 
(Chapman et al. 2004), assuming that the radio emission 
is taken to trace the starburst. Blain et al. (2004b) and 
Chapman et al. (2004) stress the need to understand the 
SEDs of these systems, specifically, in regards to ascer- 
taining whether a population of hot SMGs, which would 
be undetected in current sub-mm surveys, would con- 
tribute significantly to the infrared emission of z ~ 2 
galaxies. The most direct observational probe of the 
rest-frame far-infrared of high redshift sub-mm selected 
galaxies is the set of SHARC-2 350 /im observations of 
z ~ 2 systems obtained by Kovacs et al. (2006). CO 
observations indicate that these galaxies have large gas 
reservoirs (M gas ~ 10 10 - 10 n M Q ) (Greve et al. 2004, 
Neri et al. 2003). SMGs have also been studied through 
high resolution (~ 1") millimeter imaging and CO obser- 
vations (Genzel et al. 2003, Tacconi et al. 2006) to reveal 
large gas masses; Tacconi (et al. 2006) infer from their 
observations that their sample of SMGs do not have ex- 
tended starbursts and suggest that SMGs are likely to be 
scaled-up, more gas-rich versions of local Ultra Luminous 
Infrared Galaxies (ULIRGs), the dusty infrared-bright 
galaxies with Lg ^,m-iooo /jm ^ 1 X lO 12 ^© discovered 
in large numbers by IRAS (Soifer et al. 1984; 1987). 
SMGs have been proposed as candidates for the progen- 
itors of the most massive spheroids in the local universe 
(Lilly et al. 1999). There is also some indication that 
SMGs are a clustered population (Blain et al. 2004a), 
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with similar correlation lengths for SMGs and quasars 
(Croom ct al. 2005). 

These recent intriguing set of observations have 
prompted the development of various models to fit the 
spectral energy distributions (SEDs) of SMGs, and to 
infer an evolutionary scheme for this important class of 
galaxies. On the interpretative end, Farrah et al. (2002) 
suggest, on the basis of axially symmetric radiative trans- 
fer calculations which incorporate a model of evolving 
HII regions, that the starburst dominates in its con- 
tribution to wavelengths longer than rest-frame far-IR. 
Rowan-Robinson (2000) argues that high-redshift Hyper 
Luminous Infrared Galaxies (HLIRGs) have star forma- 
tion rates greater than 1000 M Q /yr, and are candidates 
for primeval galaxies undergoing a burst of star forma- 
tion, rather than merging systems. Efstathiou & Rowan- 
Robinson (2003) find that "cirrus" axisymmetric models, 
i.e., extended emission from dusty envelopes with low op- 
tical depths, better describe SCUBA sources, and that 
they may be more akin to optically selected high red- 
shift galaxies than obscured starbursting galaxies. Baugh 
et al. (2005) have been able to reproduce the observed 
galaxy number counts at 850 fim, assuming a top-heavy 
IMF within the context of a semi-analytic prescription 
for galaxy evolution, coupled with the axisymmetric dust 
radiative transfer code developed by Silva et al. (1998). 

High-resolution imaging of local ULIRGs has revealed 
the complex morphologies of these systems (Soifer et al. 
2000, Goldader et al. 2002, Scoville et al. 2000), indicat- 
ing a merger-driven origin. This scenario is supported by 
simulations demonstrating that tidal interactions during 
a major merger cause gas inflows by gravitational torques 
(e.g. Barnes & Hcrnquist 1991, 1996), leading to nuclear 
starbursts (e.g., Mihos & Hcrnquist 1994, 1996). A re- 
cent analysis of spectroscopic data in conjunction with 
high-resolution infrared imaging (Dasyra et al. 2006) 
indicates that the majority of ULIRGs are formed in 
nearly equal mass major mergers. Unless high redshift 
ULIRGs are more quiescent than the local systems, it 
is unlikely that axisymmetric models will prove to be 
representative of SMGs. The complex dynamical behav- 
ior of merging galaxies has been successfully modeled in 
numerical simulations incorporating feedback from cen- 
tral black holes (Springcl et al. 2005a, b; Di Matteo et 
al. 2005), and shown to reproduce the relation between 
mergers and quasar populations in optical and x-ray ob- 
servations (Hopkins et al. 2006b, c). Chakrabarti ct 
al. (2006a) employed these simulations of major merg- 
ers with black hole and starburst-driven feedback (Cox 
et al. 2006a) to calculate the infrared emission of local 
ULIRGs and LIRGs using a self-consistent fully three- 
dimensional radiative transfer code. These calculations 
reproduced observed trends such as the empirical warm- 
cold IRAS classification (de Grijp et al. 1985), wherein 
energetically active AGN are found to be correlated with 
high F 25 fj,m/F 60 //m colors (F 25 i_im/F 60 fim > 0.2), 
and demonstrated that these trends are directly driven 
by feedback processes. 

We build upon these recent developments here by cal- 
culating panchromatic SEDs, spanning optical to mil- 
limeter wavelengths, of simulations of gas-rich major 
mergers by using a self-consistent three-dimensional ra- 
diative transfer code, which treats the scattering, ab- 
sorption, and reemission of photons from dust grains 



(Chakrabarti & Whitney, 2006, in preparation). Stellar 
spectra are specified by Starburst 99 calculations (Lei- 
therer et al. 1999, Vazquez & Leitherer 2005). We con- 
duct a multiwavelength analysis of correlations from rest- 
frame near- infrared to sub-mm bands, as well as infrared 
X-ray correlations. We also compare our simulated SEDs 
with recent far-IR observations of SMGs by Kovacs et al. 
(2006). Our goal in this paper is to highlight the physical 
and dynamical basis for the multiwavelength correlations 
that we predict. In particular, our multiwavelength dy- 
namical approach allows us to cast trends in color-color 
space directly in terms of color evolution as a function 
of time, or of color evolution as a function of the relative 
luminosities from the black hole and the stars. While we 
investigate the emergent SEDs of SMGs across a wide 
wavelength range in this paper by analyzing a number 
of simulations (and different realizations of these simu- 
lations by varying the IMF and dust opacity curve, as 
we describe in §2), we defer a derivation of the infrared 
luminosity function to a future paper, as this depends 
on a full cosmological model of merger rates, and is out- 
side the scope of this paper (Chakrabarti et al. 2006c, in 
preparation) . 

The organization of this paper is as follows. In §2, we 
review the merger simulations and the translation from 
the Smooth Particle Hydrodynamics (SPH) information 
to the spatial grid that we use to do the radiative trans- 
fer calculations. In §3, we review our radiative transfer 
methodology and the dust model we have used, along 
with our adopted model for PAH emission. §4 presents 
our results, beginning in §4.1 which gives the star forma- 
tion rates, the evolution of the 850 /im fluxes, and the 
Mbh — M stal relation for SMGs. In §4.2, we present the 
simulated IRAC color-color plot in the rest-frame, and 
explain the clustering in this plot. In §4.3, we make a 
number of predictions infrared X-ray correlations, which 
can be tested empirically by obtaining a large sample of 
observations in a narrow wavelength range. We present 
photo albums spanning the lifetime of SMGs during the 
active phase in §4.4. In §5, we compare our SEDs from 
the simulations to observed data, and present the IRAC 
color-color plot in the observed frame for galaxies at 
z = 2, and show a plot depicting the variation of the 
850 fim fluxes as a function of X-ray luminosity. We in- 
troduce a simple classification scheme for SMGs on the 
basis of the Lir/£ x ratios which also corresponds to an 
evolutionary scheme. We conclude in §6. 

2. MERGER SIMULATIONS & SPECIFICATION OF 
INTERSTELLAR MEDIUM 

We employ a new version of the parallel TrceSPH 
code GADGET- 2 (Springel 2005), which uses an entropy- 
conserving formulation of smoothed particle hydrody- 
namics (Springcl & Hcrnquist 2002), and includes a sub- 
resolution, multiphase model of the dense interstellar 
medium (ISM) to describe star formation (Springel & 
Hernquist 2003a, henceforth SH03; see also Springel & 
Hernquist 2003b) . Black holes are represented by "sink" 
particles that accrete gas, with an accretion rate esti- 
mated using a Bondi-Hoyle-Lyttleton parameterization, 
with an upper limit equal to the Eddington rate (Springel 
et al. 2005b). The bolomctric luminosity of the black 
hole is then Lboi = e Y Mc 2 , where e r = 0.1 is the radia- 
tive efficiency. We further allow a small fraction (~ 5%) 
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TABLE 1 

SMG Simulations 



Simulation Progenitor Redshift Orbital Orientation V T ; r (km/s) 



cl60(vc3e) e 160 

hl60(A3) h 160 

h226(A4) h 226 

h320(A5) h 320 

h500(A6) h 500 

c226(A4e) e 226 

c320(A5e) e 320 

c500(A6e) e 500 

z3c270(z3A4.5e) 3 h 270 

z3h500(z3A6) 3 h 500 



of Lboi to couple dynamically to the gas as thermal en- 
ergy. This fraction is a free parameter, determined in Di 
Matteo et al. (2005) by matching the Mbh-c relation. 
We do not attempt to resolve the gas distribution imme- 
diately around the black hole, but instead assume that 
the time-averaged accretion can be estimated from the 
gas on the scale of our spatial resolution, which is a few 
tens of parsecs in the best cases. 

We adopt the ISM prescription developed by 
Chakrabarti et al. (2006a) in their study of the infrared 
properties of local ULIRGs. Owing to the limited res- 
olution of current numerical simulations, the interstel- 
lar medium (ISM) must be treated in a volume-averaged 
manner. In particular, the simulations discussed here 
adopt the multiphase model of SH03, in which individ- 
ual SPH particles represent a region of the ISM that con- 
tains cold clouds embedded in a diffuse hot medium. Be- 
cause the SPH calculation uses only the volume-averaged 
pressure, temperature, and density to evolve the hydro- 
dynamics, this model does not provide specific informa- 
tion regarding the cold clumps. Since these cold clumps 
harbor the dust that produces the infrared emission, we 
must adopt a model which determines their number, den- 
sity, size, and location. 

As the first step to specify the properties of the cold, 
dense gas, we assume that each SPH particle contains 
one dense molecular cloud at its center. Because we know 
the temperature and density of the volume-average SPH 
particle, assuming a temperature (say, 10 3 K, as was as- 
sumed in SH03) for this cold cloud and that it exists 
in pressure equilibrium with the hot phase we readily 
attain its density. However, observations of molecular 
clouds indicate that their pressure is dominated by tur- 
bulent motions (Blitz et al. 2006; Solomon et al. 1997), 
rather than thermal energy, and thus neglecting this fea- 
ture yields very dense clouds, with small volume filling 
factors. In the work presented here, we take the tur- 
bulent pressure to be ~ 100 times that of the thermal 
pressure; the turbulent velocity as estimated from line 
widths is of order 10 times the thermal sound speed in 
star forming regions, which would lead to a factor of 
~ 100 for the ratio of the turbulent to thermal pressure 
(Plume et al. 1997). With this assumption, the density 
is much lower than without the turbulent support and its 
volume filling factor is greatly increased. (Our approach 
is similar to that employed by Narayanan et al. [2006] 
in their study of the evolution of the molecular gas in 



mergers.) 

We have considered a diverse range of simulations in 
this paper, of varying mass, orbital inclination, and pro- 
genitor redshift. We give a summary of the simulations 
analyzed in this paper in Table 1, including the pro- 
genitor redshift of the two galaxies, the orbital orien- 
tation, which is denoted "h" for a co-planar merger, and 
"c" for an arbitrary inclination, and the virial velocity. 
The names that we use to refer to the simulations along 
with names used in previous papers (i.e., Hopkins et al. 
2005b) given in parentheses. A detailed description of 
these models is given in Cox et al. (2006a). The orbital 
orientation denoted "e" corresponds to 9\ = 30, 4>i = 60 
for the first disk, and O2 = —30, <j>2 = 45 for the second 
disk (also in Table 1 of Cox et al. 2006a). We construct 
the disk progenitors at redshift z in a manner similar to 
Robertson et al. (2006) which follows the Mo, Mao & 
White (1998) formalism in determining the scale lengths 
of the disks, and the disk concentrations are scaled fol- 
lowing Bullock et al. (2001). 

3. RADIATIVE TRANSFER METHODOLOGY 

We use a self-consistent three-dimensional Monte Carlo 
radiative equilibrium code (Chakrabarti & Whitney, 
2006, in prep) to calculate the emergent SEDs and images 
from the merger simulations as a function of evolution- 
ary state. We review some of the main points here. This 
code incorporates the Monte Carlo radiative equilibrium 
routine developed by Bjorkman & Wood (2001) to solve 
for the equilibrium temperature of the dust grains, in a 
manner similar to that implemented by Whitney et al. 
(2003). Individual photons are tracked directly in the 
Monte Carlo scheme. When a photon is absorbed by a 
grid cell, it may raise the temperature of that grid cell so 
that it emits infrared radiation; the emissivity of this grid 
cell depends on the temperature. The Bjorkman & Wood 
(2001) algorithm corrects the temperature in a grid cell 
by sampling the new photon frequency from the differ- 
ence of two emissivity functions - that from the previous 
temperature and that computed with the new tempera- 
ture. This is not an explicitly iterative scheme; a suffi- 
ciently large number of photons must be absorbed by the 
grid cells so that they eventually relax to an equilibrium 
temperature. Since the density within a grid cell is con- 
stant, one can easily integrate the optical depth through 
the grid to account for the attenuation of the photons. 
This code also treats scattering, representing the scatter- 
ing phase function as a Henyey-Greenstein phase func- 
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tion (see Whitney et al. 2003; Whitney & Wolff 2002 
for a detailed description of this). Since the scattering 
efficiency is proportional to 1/A 4 , it is particularly rele- 
vant for the shorter wavelengths. Whitney et al. (2003) 
carried out convergence studies to compare their results 
to a set of benchmark calculations developed for spher- 
ically symmetric codes. Chakrabarti & Whitney (2006, 
in prep) have carried out both photon and grid resolu- 
tion convergence studies to verify that photon numbers 
of order 10 8 yield converged results from mm to near-IR 
wavelengths for the three-dimensional grids we have used 
here, which follow the clump size distribution (which re- 
sult from our ISM prescription discussed in §2) and will 
at a minimum, resolve (i.e., by ~ 2 cells) the median 
of the clump size distribution; the grid resolution was 
varied by an order of magnitude to find converged re- 
sults. For fine grid structures, it is necessary to increase 
the photon number proportionately to ensure that each 
grid cell receives enough photons to reach an equilibrium 
temperature (i.e., the photon and grid resolution stud- 
ies are not independent). Off- nuclear sources of radia- 
tion are allocated photons in proportion to the fraction 
of bolometric luminosity they contribute. The radiative 
equilibrium temperature calculation needs to be refer- 
enced to the total bolometric luminosity of all of the 
sources of radiation. For large three-dimensional grids, 
it is crucial to allocate tens of millions of photons for 
the grid cells to reach an equilibrium temperature; oth- 
erwise, the emissivity of the dust may be artificially low- 
ered. This code incorporates a "peeling-off" algorithm 
(Yusef-Zadeh, Morris & White 1984; Wood & Reynolds 
1999) to calculate the images, which are then convolved 
with broadband filter functions. 

We model the intrinsic AGN continuum spectrum fol- 
lowing Hopkins, Richards & Hernquist (2006) (HRH), 
which is based on optical through hard X-ray observa- 
tions (Elvis et al. 1994, George et al. 1998, Perola et 
al. 2002, Telfer et al. 2002, Ueda et al. 2003, Vignali 
et al. 2003), with a reflection component generated by 
the PEXRAV model (Magdziarz & Zdziarski 1995). The 
HRH spectrum is similar to that developed by Marconi 
et al. (2004), but it is more representative of the shapes 
of typical observed quasars; it includes a template for 
the observed "hot dust" component longward of 1 /im. 
Since we calculate the dust emission self-consistently, we 
enforce a Rayleigh- Jeans cutoff beyond 1 /im in the in- 
put AGN spectrum. In our models here, the presence of 
dust emission at wavelengths longer 1 fim is derived from 
our radiative transfer calculations. We have performed 
Starburst 99 calculations (Leitherer et al. 1999, Vazquez 
& Leitherer 2005) to calculate the stellar spectrum and 
bolometric stellar luminosity, taking as input the age, 
mass, and metallicity of the stars from the simulations, 
for both the Kroupa and Salpeter IMFs. 

We use the Weingartner & Draine (2001) (henceforth 
WD01) Rv = 5.5 dust opacity model, both with and 
without the addition of ice mantles. The inclusion of ice 
mantles (which lead to an effective increase in the opacity 
of ~ 2 for T < 100 K above which ices sublimate, e.g., 
Pollack et al. 1994) is motivated by spectroscopic studies 
that have identified ice spectral features in protostcllar 
environments and in dusty galaxies (Allamandola et al. 
1992; Spoon et al. 2002). The mass fraction of dust is 
equal to 1/105.1 for solar abundances (WD01). WDOl's 



models have been shown to reproduce the observed ex- 
tinction curves for the Milky Way, as well as regions of 
low metallicity, such as the LMC and the SMC. The opac- 
ity normalization per gram of dust, n\ 01 is equal to 0.275 
for A = 100 ^im, where 5 is the dust-to-gas ratio rela- 
tive to solar, which we take to be unity. Dunne & Eales 
(2001) and Klaas et al. (2001) found' that the dust-to- 
gas ratio for a large sample of ULIRGs is comparable to 
Milky Way values, when they fit two-temperature black- 
bodies to the far-IR SEDs. Wilson et al. (in preparation) 
also find similar results when they use multi-temperature 
blackbodics to fit the SEDs. Previous work, based on 
fitting single temperature blackbodies had found lower 
dust-to-gas ratios (Dunne et al. 2000) by a factor of 
two. Further, this dust model has been used successfully 
in fitting the observed SEDs of ULIRGs and SMGs by 
Chakrabarti & McKee (2005) and Chakrabarti & McKee 
(2006, in preparation). We adopt a phenomenological 
model of PAH emission here. If there is a UV photon 
of sufficient magnitude to excite PAH emission, then we 
add on an observed PAH template (adopted from Do- 
pita et al. 2005) to the SED, where the PAH strength is 
scaled to the 200 ^m flux as in Dale & Helou (2002). 

4. RESULTS & ANALYSIS 

We plot in Figure 1 the star formation rates for the 
simulations analyzed in this paper. SH03 implement an 
algorithm for star formation that is similar to the Ken- 
nicutt (1998) relation, i.e., Ssfr oc Sg a t- For progen- 
itors at a fixed redshift and gas fraction, the star for- 
mation rate is proportional to (essentially a linear power 
of) the virial velocity, since M vir = V^ 3 ir /(10 GH{z)) and 
i? vir = Kir/(10 GH(z)) which gives S S fr oc for 
progenitors at the same redshift, with fixed gas fraction. 
The sequence of plots in Figure 1 (a) and (b) shows that 
the peak star formation rate scales with virial velocity 
for progenitors at a given redshift. Disk progenitors at 
redshift z are constructed in a manner similar to that 
described by Robertson et al. (2006); the size of the 
disks are set by the Mo, Mao & White (1998) formalism 
wherein the disk contains a fraction of the total galaxy 
angular momentum equal to its mass fraction. The disk 
scale length is then determined by the galaxy spin (which 
is set to a constant that is motivated by results from cos- 
mological N-body simulations) and the dark matter halo 
concentration, C v ; r , which scales essentially as 1/(1 + z) 
(with a weak dependence on the virial mass, e.g., equa- 
tion 3 in Robertson et al. 2006). This dependence of the 
concentration on the redshift offsets the 1/H(z) depen- 
dence of i?vir (i.e., when substituting this dependence of 
the concentration in eq. 28 of Mo et al. 1998) so that 
the disk scale radius decreases by about 20 % to z = 3 
while there is a factor of ~ 5 loss in mass. Therefore, the 
star formation will depend primarily on the mass of the 
disk (and will decrease proportionately) for progenitors 
at redshift z. The scaled z = 3 simulations, for a given 
virial velocity (e.g. z3h500 (z = 3) in Figure lc compared 
to h500 in Figure la) have lower star formation rates, as 
they are less massive. We also discuss later in §4.4 that 
the simulations with progenitors at z = 3 have more com- 
pact morphologies than the simulations with disk prop- 
erties of z = simulations. Inferred star formation rates 
of SMGs vary from the extremely high star formation 
rates of cirrus models, SFR > 1000M Q /yr (Efstathiou & 




Fig. 1. — Star Formation Rate (in M@/yr) vs. time for the merger simulations analyzed in this paper. Left panel shows the z = 
co-planar mergers of varying virial velocity; the number refers to the virial velocity of the simulations in km/s; middle panel shows the 
2 = non co-planar mergers; right panel shows the scaled z = 3 mergers ("e" denotes arbitrary inclination and "h" is co-planar) 



Rowan- Robinson 2003; Rowan- Robinson 2000), to about 
5OOM0/yr if a factor of two is used to account for the 
AGN contribution to the infrared luminosity, along with 
Kennicutt's (1998) relation for the star formation rate 
and the infrared luminosity (Genzel et al. 2003); infer- 
ring star formation rates without including the contribu- 
tion of the AGN leads to star formation of the order of 
several thousand solar mass per year (Ivison et al. 1998). 
In our simulations, during the active phase, the star for- 
mation rates are on average ~ 500 — 1000Af©/yr and 
produce bright SMGs as we discuss below. 

Henceforth, we will illustrate specific features of the 
merger simulations for representative samples of simula- 
tions. We show in Figure 2 the observed 850 /im flux as 
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Fig. 2. — i*850 /im (in observed frame for galaxies at z = 2, ss 
15.5 Gpc) vs. time for the SMGs of varying virial velocity (which 
corresponds to varying galaxy mass for galaxies with progenitors at 
the same redshift) in our sample; solid line is the el60 simulation, 
dotted line is e226, dashed line is e320. Solid horizontal line marks 
^850 fim > 1 mjy, the empirical definition of SMGs. These all 
have the same orbital inclination. 



a function of time for el60, e226, and e320, which span 
the range from virial velocities of 160 km/s to 320 km/s. 
The 850 /jm flux is shown for galaxies at a luminosity 
distance of L>l ~ 15.5 Gpc (this corresponds to the rest- 
frame 283 /im flux at z = 2). SMGs are an empirically 
defined class of galaxies - based on the observed 850 /im 
flux; systems with observed Fg5o /jm ^ 1 mJy are desig- 
nated as sub-mm bright (Ivison et al. 2000; Blain et al. 
1998). The solid horizontal line demarcates this criterion 
(■^850 /im ^ 1 mJy) on our flux versus time plots. (This 
empirical definition of sub-mm bright may well change 
owing to improved sensitivity of future sub-mm surveys, 
i.e., SCUBA-2) As is clear, the peak 850 /im flux in- 
creases as a function of virial velocity (which corresponds 
to galaxy mass for galaxies with progenitors at the same 
redshift). The more massive galaxies show a more rapid 
and pronounced decline in their sub-mm flux as a func- 
tion of time as these systems are losing gas content both 
owing to AGN feedback clearing out the dust and gas in 
a more violent fashion (Hopkins et al. 2006d, Hopkins & 
Hernquist 2006), as well as to forming stars at a faster 
rate. 

Figure 3 shows the observed 850 /im flux as a func- 
tion of time for galaxies with varying orbital inclination 
for e226, h226, e320, and h320. While there is no clear 
difference between the peak fluxes of the galaxies shown 
here of varying orbital inclination, the co-planar mergers 
(shown here as the solid lines) do have a sharper de- 
cline in their 850 /im flux as a function of time. Finally, 
Figure 4 contrasts the time evolution of the 850 /im flux 
when the WD01 opacity (solid line) is used instead of the 
WD01 model with inclusion of ice mantles (dotted line). 
Since the observed 850 /im flux is nearly optically thin, 
the effect of adding ice mantles corresponds to increasing 
the 850 /im flux. Figure 5 shows the e226 and e320 simu- 
lations with the Salpeter (solid line) and Kroupa (dashed 
line) IMFs. Using the Kroupa IMF leads to a slight in- 
crease in the observed 850 /im flux (by about ~ 30%). 
Contrasting Figure 1 with Figures 2-5 demonstrates that 
the observed 850 /im flux does track the star formation 
rate quite closely. We return to this point in a forthcom- 
ing paper where we analyze a larger set of simulations 
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and give fitting relations between the 850 /im flux and 
relevant parameters, such as the star formation rate. 

Figure 6 shows the time evolution of the observed 
850 fim fluxes for the el60 and hf 60 simulations - which 
are designed to have virial velocities representative of 
local ULIRGs (the el60 simulation has 40% gas frac- 
tion initially, while the hl60 simulation has 100 % at the 
start of the simulation). These simulations have sub- mm 
fluxes similar to those of local LIRGs and ULIRGs. Av- 
eraging the 850 fim fluxes in Table 3 of Dunne & Eales 
(2001) gives 183 mJy, with Arp 220 topping the list at 
832 mJy. Such objects at Dl ~ 77 Mpc would be bright 
in the sub-mm for a longer fraction of their lifetime than 
at z = 2, as has been shown previously in this section. 
We discuss in Chakrabarti et al. (2006c, in preparation) 
that the time variation of the 850 /im fluxes leads to 
a simple interpretation for the observed 850 fim num- 
ber counts, as well as the clustering properties of SMGs, 
similar to the characterization of the quasar luminosity 
function by Hopkins et al. (2005a-c, 2006a), and the 
weak luminosity dependence of clustering properties of 
quasars, as discussed by Lidz et al. (2006). Although we 
do not derive the sub-mm luminosity functions in this pa- 
per, it is clear that simulations of gas-rich major mergers 
naturally lead to the production of SMGs. 

We show in Figure 7 the evolution of MBn/M sta r as 
a function of time for e226 (red line), e320 (blue line) 
and e500 (green line). At late times, we see that these 
simulations tend to the observed M sta r — Mbh relation 
for spheroids, i.e., the mass of the black hole is ~ 0.1% 
of the spheroid mass (Magorrian et al. 1998). Borys et 
al. (2005) discuss the relationship between the stellar 
and black hole mass for SMGs at a median redshift of 
2. Their X-ray observations indicate that about half of 
their sample of 13 SMGs host AGN. They derive the re- 
lation between the black hole mass, from the observed 
X-ray luminosity under the assumption of Eddington- 
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Fig. 3. — i*850 (im (in observed frame for galaxies at z = 2, Dl ~ 
15.5 Gpc) vs. time for the SMGs of varying orbital inclination for 
h226, e226, h320, and e320. The solid line shows the co-planar 
cases: h226 and h320, and the dotted line the non co-planar cases 
c226 and c320. 
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Fig. 4. — ^850 fim (in observed frame for galaxies at z = 2, 
Dl ~ 15.5 Gpc) vs. time for the SMGs of varying dust opacity 
shown here for e226, with and without the inclusion of ice mantles. 
Solid line has the WD01 opacity, while the dotted line includes ice 
mantles. 
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Fig. 5. — ^850 urn (in observed frame for galaxies at z = 2, 
Dl ~ 15.5 Gpc) vs. time for for e226 and e320, shown here for 
Salpctcr (solid line) and Kroupa (dashed line) IMFs. 



limited accretion, and stellar mass, from the rest-frame 
K band observations, to find that the SMGs fall about 
two orders of magnitude below the local relation. The 
SMG phase for the simulations shown in Figure 7 oc- 
curs around t ~ 0.5 h _1 Gyr - such systems would indeed 
lie below the local relation, but not by two orders of 
magnitude. We show the change in MBu/M s t ar when 
we obtain the black hole masses assuming Eddington- 
limitcd growth (the dotted lines in this figure) - which 
would then lead to the appearance of such objects falling 
about two orders of magnitude below the local relation. 
While we find that the black hole masses would grow 
perhaps by a factor of ~ 5 between the SMG phase at 
z = 2 and the present-day, it is possible that Borys et 
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Fig. 6. — i*850 fim vs. time for galaxies designed to be rep- 
resentative of the local population (el60 and hl60), placed at 
D L = 77 Mpc. 



al.'s (2005) assumption of Eddington-limited growth has 
led to an overestimate of this factor by more than an 
order of magnitude. 

4.1. IRAC Correlations 

We plot in Figure 8 an IRAC color-color plot in the 
rest-frame. This IRAC color-color plot was proposed by 
Lacy et al. (2004) to identify AGN; the AGN are the 
sources that fall within the dashed lines. We recover 
the general trends of the observed IRAC color-color plot 
from Spitzer Space Telescope's First Look Survey (FLS) 
- namely, the clustering of the fluxes in the lower left 



hand corner of the plot (both colors being bluer), as well 
as a small fraction of these sources falling on the redder 
sequence. We show a much smaller sample here than 
that shown by Lacy et al. (2004), which depicts tens of 
thousands of sources, while we consider ten simulations, 
at about 20 different snapshots during the active phase 
on average. Specifically, of the 16,000 objects shown in 
the Lacy plot, 2000 are likely to contain AGN, while 
20 were robustly identified to be obscured AGN; about 
eight out of about of about two hundred points fall in the 
dashed region from our simulated data. (For each time 
snapshot in the simulation, the median flux over viewing 
angles constitutes a single data point). 

The power of the dynamical approach afforded by cal- 
culating fluxes from the simulations is that we can un- 
fold this color-color plot to a color vs. time plot for 
each axis, as shown in the following Figure 9. As this 
plot shows, the clustering that is present in the observed 
color-color plot, which we reproduce also in the simu- 
lation IRAC color-color plot, is naturally explained in 
our models by the SMG spending more of its lifetime 
in that region of color space. The low F% um./Fi.5 fim 
and F5.8 fim/Fs.e colors correspond to the Rayleigh- 
Jeans tail of attenuated starlight (modulated by the 
albedo) influencing the IRAC colors. (It is a simple exer- 
cise to show that the Rayleigh- Jeans tail of stellar spectra 
leads to log(F 8 / im/F 4 .5 /urn) and log(F 5 . 8 ^ m /F 3 .6 jum) 
roughly of order -0.5, independent of the details of the 
stellar spectrum.) Figure 10 plots the IRAC colors as 
a function of LBH/^stars - this shows that if the IRAC 
colors fall within the dashed region demarcated by Lacy 
et al. (2004), the source is certainly an energetically ac- 
tive AGN; i.e., the bolometric luminosity from the black 
hole is greater than or comparable to the luminosity 
from the stars. However, there are also cases, when the 
AGN's luminosity is greater than or comparable to the 
stars, and yet it is not in the AGN-demarcated region of 
color-color space. For instance, most of the points which 
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Fig. 7. — A^Bii/A^star versus time for e226 (red line), e320 (shown 
in blue), and e500 (shown in green). The dotted lines correspond to 
deriving the black hole mass for Eddington limited accretion. This 
shows that deriving black hole masses under the assumption of 
Eddington-limited accretion underestimates the black hole masses 
by a factor of ~ 5 during the SMG phase (t ~ 0.5h _1 Gyr). 



Fig. 8. — Correlations in rest-frame IRAC bands, shown for all 
the simulations analyzed in this paper. Dashed lines demarcate 
color criteria to pick out obscured AGN sample by Lacy et al. 
(2004). 
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Fig. 9. — Rcst-framc IRAC Color-Color plots as a function of time. The clustering in the lower left hand corner of the IRAC color-color 
plot owes to the simulated galaxies spending most of their time in that phase. 



have blue colors in F 5 . 8 ^m/-F3.6 /im and red colors in 
^8 j um/-F 1 4.5 fj,m do have Lbh > istars, as well as a few 
of the points that have blue F 5 . 8 ^m and blue 

^8 j um/-F 1 4.5 fim- (For the most part, though, sources 
that have blue colors for both axes are dominated by the 
stellar luminosity.) This variation seems to owe to the 
effects of larger and more extended starbursts in some 
of the co-planar mergers. Therefore, while we confirm 
Lacy's et al. (2004) IRAC correlations in general, our 
interpretation of the blue-red plume, as well as sources 
that lie close to the demarcating line for AGN, is differ- 
ent from their analysis. We suggest that some of these 
sources may also have energetically active AGN, which 
could be detected in deep X-ray surveys. 

We emphasize that this interpretation of the IRAC 
color-color plot is most appropriate when viewed in the 
rest-frame, as the clustering properties and relative frac- 
tion in the AGN demarcated region can then be directly 
correlated with dynamical properties of these galaxies, 
such as the relative amount of time the system spends in 
some region of color-color space and the relative contri- 
bution from the black hole luminosity. One can interpret 
the FLS IRAC color-color plot roughly in the manner 
outlined above as well since Lacy et al. (2004) note that 
the median photometric redshift of their candidate AGN 
is ~ 0.3. We discuss later in §4.5 that the continuum is 
not as sensitive to redshifting the SED as the PAH fea- 
tures - hence the continuum in the observed frame is quite 
similar (for a range of redshifts at least, as we discuss in 
§4.5) to the rest-frame color-color plot. Nonetheless, the 
dynamical interpretation of trends in color-color space is 
most easily understood in terms of rest-frame quantities. 

The recent, comprehensive work by Barmby et al. 
(2006) corroborates our suggestion regarding the blue- 
red plume containing energetically active AGN. Specifi- 
cally, Barmby et al. (2006) find that the infrared proper- 
ties of X-ray sources in the Extended Groth Strip are very 
diverse - about 40% of the X-ray detected sources have 
red power-law SEDs in the 3.6 — 8.0 ^m IRAC bands, 
while another 40% have blue SEDs. Thus, while the 
presence of sources in the AGN demarcated region in 
the IRAC color-color plot is a sure sign of an energeti- 
cally active AGN, the converse is not necessarily true (as 



Figure 10 has shown). The diversity of infrared proper- 
ties of X-ray detected sources is further highlighted by 
Rigby et al. (2004), who found that X-ray hard AGNs 
are not preferentially infrared-brighter; Alonso-Herroro 
et al. (2004) found that sources in the Lockman Hole 
with similar selection criteria have a variety of optical/IR 
spectral types. Therefore, while AGN identification us- 
ing the IRAC color-color selection is generally effective, it 
may undercount the number of energetically active AGN 
as it would miss the bluer region. 

We show in Figure 11 our simulated IRAC color-color 
plot with inclusion of our phcnomcnological model of 
PAH emission. The PAH region, as shown in this fig- 
ure, is qualitatively in agreement with Sajina et al.'s 
(2005) results, who use a combined template of stellar 
emission, a power-law for the mid-IR spectrum, and a 
PAH template to simulate an IRAC color-color plot. In- 
terestingly, our color-color plot for the continuum (as 
shown in Figure 8) - which is derived from self-consistent 
three-dimensional radiative transfer solutions - is in bet- 
ter agreement with the observed FLS IRAC color-color 
plot shown in Lacy et al. (2004) than Sajina et al.'s 
(2005) results. Lacy et al. (2004) note that the median 
photometric redshift of the candidate AGN is ~ 0.3 in 
the FLS color-color plot. The location of the PAH re- 
gion in the color-color plot is a sensitive function of red- 
shift. While in the rest-frame the PAH region extends 
into (and bends towards) the AGN demarcated region 
(as we have shown here), redshifting the PAH template 
to z > 0.1 causes the PAH region to bend into the bluer 
region - resulting in the appearance of the "bunny-ear" 
shape seen in FLS IRAC color-color plot, as has been 
noted by Sajina et al. (2005). We show later in §4.5 an 
IRAC color-color plot in the observed frame for galax- 
ies at a range of redshifts - from z w 0.3 — 2, which do 
indeed result in the characteristic "bunny-ear" shape for 
the z = 0.3 sample. 

We note that that the scaled z = 3 simulations do not 
occupy a region of color-color space (in the rest-frame) 
that is distinct from that of the z = simulations. This 
suggests that at least qualitatively, the appearance of (in 
particular, the clustering of the colors) simulated galaxies 
on color-color plots is a generic feature of simulations of 
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Fig. 10. — Rest-frame IRAC Color-Color plots as a function of Z/Bll/Lstars. The clustering in the lower left hand corner of the IRAC 
color-color plot is positively correlated with the stars dominating in their contribution to the total bolomctric luminosity. The redder colors 
are positively correlated with the black hole being energetically active, i.e., Lbh ^ L s tars- There are however cases when the colors are 
blue and Lbh ^ L s tars (see discussion in §4.1). 



1.5 
1.0 
0.5 



g 1 0.0 



-0.5 



1.0 




1.0 -0.5 0.0 0.5 1.0 1.5 

l0 g( L 5. 8 / L 3.6) 



Fig. 11. — Correlations in rcst-framc IRAC bands, with inclusion 
of PAH emission model (squares designate the combined spectrum, 
including the continuum and PAH model). Dashed lines demarcate 
color criteria to pick out obscured AGN sample by Lacy et al. 
(2004). 



major mergers where black hole feedback has a significant 
and sharp effect on the evolution of the system, but which 
evolve otherwise passively. Extrapolating this trend to 
high redshift systems suggests that simulations of high- 
rcdshift systems (z > 6) would show a similar behavior 
in color-color space. 

4.2. Infrared X-ray Correlations 

A large population of obscured AGN is discerned from 
the hard X-ray background (Comastri et al. 1995); these 
systems would re-radiate a significant fraction of their 
luminosity into the infrared, in principle explaining the 
spectral shape of the X-ray background (see, e.g., Hop- 
kins et al. 2006a, HRH). Therefore, infrared X-ray cor- 
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Fig. 12. — (a) Top panel - predictions for infrared luminosity 
(integrated from 8 (im to 1000 (im) vs. the hard X-ray, (b) Middle 
panel - predictions for 70 /im rest-frame luminosity vs. the hard 
X-ray correlations, (c) Bottom panel - predictions for 24 ^m rest- 
frame luminosity vs. the hard X-ray luminosity. 



relations hold promise for being particularly useful di- 
agnostics for discerning and understanding the nature 
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of both energetically active AGN, and dormant AGN 
approaching their accretion phase in merging galaxies. 
We show in Figure 12 (a) the infrared luminosity ver- 
sus the hard X-ray luminosity (from 2-10 keV) from our 
set of simulations. In our model, the hard X-ray lumi- 
nosity derives from the black hole luminosity which we 
compute following the model of Marconi et al. (2004). 
On the high X-ray luminosity end (L x > 10 9 £q), the 
contribution from the black hole to the X-ray luminos- 
ity is significantly higher than that from star forma- 
tion; even moderately high values of the star formation 
rate (~ 100 M /yr) lead to x-ray luminosities of order 
~ 10 8 £©/yr (Gilfanov et al. 2004). X-ray emission 
from hot gas is also negligible on the high x-ray luminos- 
ity end - from analyzing the x-ray emission from hot gas 
in merger simulations, Cox et al. (2006b) find that the 
maximum x-ray emission would not exceed ~ 10 s Lq. 

Figure 12 (a) shows that there is a clear correlation 
between the the total infrared luminosity (which we get 
from integrating the SED from 8 /im to 1000 ^m) and 
the hard X-ray luminosity, which extends over nearly four 
orders of magnitude, with an increase in the scatter at 
low X-ray luminosities. This makes intuitive sense - in 
our model, low X-ray luminosities correspond to the pre- 
(and post-) accretion phase, before (and after) the black 
hole has become energetically active (Lbh ^ L stal )- In 
the low X-ray luminosity phase, the infrared luminosity 
is partly powered by distributed sources of radiation; i.e., 
stars, and partly by a heavily enshrouded black hole. On 
the most X-ray luminous end, the black hole dominates 
in its contribution to the total bolometric luminosity - 
this case is essentially analogous to a central source of 
luminosity powering the infrared radiation. Hence, there 
should be a correlation between the luminosity of the cen- 
tral source and the reprocessed radiation in this phase. 
These systems, at their peak infrared luminosity, usually 
exceed ULIRG, and sometimes HLIRG (L m > 10 13 i Q ) 
luminosities. On average, the X-ray luminosities are of 
order 10 10 L Q , though there are a few on the very bright 
X-ray end (L x - Iay ~ 10 n £ Q — 10 12 L©) - which corre- 
spond to the phase when the black hole dominates the 
bolometric luminosity. Alexander et al. (2005a) find 
that the X-ray luminosity of their SCUBA selected sam- 
ple is ~ lO 11 ^©, with only four in their sample exceed- 
ing 10 n £ Q . Alexander et al. (2005a) use the locally 
observed radio far-IR correlation to find that the X-ray 
luminosity is about 250 times lower than the infrared 
luminosity, whereas we find that it is typically ~ 100 
times lower. The radio far-IR correlation has consider- 
able scatter in luminous galaxies (Sadler et al. 2002), 
and its extrapolation to high redshifts is uncertain. As 
such, this discrepancy may in part owe to the fact that 
using the radio far-IR correlation underpredicts the in- 
frared luminosity since the contribution from the black 
hole has not been taken into account. Another possibil- 
ity is that (part of) their X-ray sample is Compton thick, 
and hence have higher intrinsic luminosities than have 
been accounted for. Multi-wavelength infrared observa- 
tions are needed to observationally derive the rest-frame 
infrared luminosities to test our prediction. 

We find that the 70 ^m luminosity is roughly forty 
times the X-ray luminosity, and displays the same behav- 
ior as the infrared luminosity (Figure lib). Figure 12 (c) 
shows that a similar behavior for the rest-frame 24 /im 



flux - with the brightest 24 /zm objects being strongly 
correlated with the hard X-ray luminosity. We stress 
that these correlations are predictions from our model, 
which can be directly tested by combining far-IR and 
X-ray data at comparable sensitivities for large samples 
of SMGs. This should be made possible by the upcom- 
ing Hcrschel mission. Source-stacking analysis may make 
testing this prediction feasible, at least in a statistical 
sense, by using existing MIPS data. 

We also show in Figures 13 (a-d) the rest-frame lumi- 
nosities, derived from the continuum only, in the IRAC 
bands plotted versus the hard X-ray flux. These plots 
display a more complex behavior than the far-infrared 
X-ray correlations. There is a positive correlation only 
at high X-ray luminosities, with the weakest correla- 
tions, i.e., increase in Lirac bands as a function of X- 
ray luminosity, seen for the shorter wavelengths. The 
different curves owe to more massive simulations hav- 
ing higher luminosities; actual observations will probe 
nearly a continua of masses and hence will not appear 
discrete. At the longer IR wavelengths, the dust opac- 
ity curve has a relatively simple (essentially power-law) 
form and the emission comes from the outer regions of 
the envelope where the dust is cooler. However, the 
rest-frame IRAC fluxes typically emanate from optically 
thick regions, which are often heated to high tempera- 
tures (T > 100K), with their emissivity modulated by 
(distinctly non power-law) spectral features in the dust 
opacity curve. (For a discussion of, and analytic expres- 
sions for, the characteristic emission radii of different fre- 
quency regimes in dust envelopes, see CM05.) Nonethe- 
less, on the bright end, the problem does simplify for the 
longer wavelength IRAC bands. The shortest wavelength 
is the least correlated with the X-ray luminosity - dust is 
rarely heated to temperatures high enough (T > 900 K) 
to produce 3.6 /im emission, except on the very bright 
end. As such, this band is generally probing the atten- 
uated Rayleigh-Jeans fall-off of the stellar light, partic- 
ularly at low X-ray luminosities. The fluxes in these 
bands are also influenced by scattering of photons by dust 
grains since the scattering efficiency increases sharply at 
short wavelengths and leads to an increase in the flux 
by factors of 3-5 (for optical depths characteristic of the 
simulations), especially as the effective optical depth of 
the envelope decreases (Chakrabarti & Whitney, 2006, in 
preparation). Moreover, they will also be influenced by 
PAH emission. Yan et al. (2005) find that 60% of their 
sources between redshifts of 1.8 and 2.6, have weak or 
no PAH emission, while about half have prominent sil- 
icate absorption lines. Both PAH spectral features and 
non-power law continuum features will heighten the rel- 
ative scatter for IRAC X-ray correlations, as compared 
to the far-IR X-ray correlations. Rigby et al. (2004) 
emphasize the diversity of 24 /im properties of X-ray se- 
lected AGN, even in a narrow redshift slice (z ~ 0.7). 
They do see some indication of the ratio of the 24 /im 
to X-ray luminosity increasing with X-ray hardness for 
their subsample at z ~ 0.7, but there is scatter even in 
this subsample. A detailed comparison with observations 
is beyond the scope of this paper - in a future paper, 
we use the SEDs from the simulations as templates to 
K-correct observed fluxes of sources with spectroscopic 
redshifts, using the redshift distribution from Chapman 
et al. (2005) for SMGs. 
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Fig. 13. — (top-left) The rest-frame 3.6 fim luminosity vs. the hard X-ray, (top-right) the rest-frame 4.5 /im vs. the hard X-ray, 
(bottom-left) rest-frame 5.8 (im vs. the hard X-ray luminosity, (bottom-right) rest-frame 8 (im vs. hard X-ray luminosity. Note that the 
clearest correlations are for the longer wavelength bands at higher X-ray luminosities (see discussion in §4.2). 



4.3. Photo Albums of the Lifetimes of SMGs 

We present photo albums of the time evolution of the 
surface brightness in the observed 3.6 /im band, as well 
as the 450 /im band during the lifetime of three SMGs 
from our sample. We adopt a luminosity distance of 
15.5 Gpc corresponding to z = 2 (Qm = 0.3, Q\ = 0.7, 
H a = 70) in these figures and take into account the in- 
strumental angular resolution. These figures should not 
be interpreted as what would be observed by SCUBA and 
Spitzer's IRAC since we have not convolved these images 
with realistic PSFs; doing so would not leave any visible 
structure. Our goal here is merely to use these images 
as a visual aid in contrasting the time evolution of the 
surface brightness and morphologies of three SMGs of 
differing orbital inclination and progenitor rcdshift dur- 
ing their lifetimes. 

Figure 14 shows a co-planar merger (h320) from its pre- 
merger phase (a), to close to the main feedback phase as 
the AGN starts to clear out the surrounding dust and 
gas (b), to the final remnant (c). Figure 15 (a-c) shows 
the corresponding images in the observed 450 ^m band 
at the same times. The galaxy becomes brighter in the 
IRAC bands as it is becoming progressively fainter in the 
longer wavelength SCUBA band. This simply reflects the 
time variation of the fluxes - as the AGN feedback clears 



out the obscuring material, more of the emitted energy 
(the SED) shifts to shorter wavelengths, and there is cor- 
respondingly less emission at longer wavelengths. This 
trend would be seen also in surface brightness maps of 
the CO emission which probe the cooler regions of the 
envelope. Another point of note is that this co-planar 
merger has a disk-like morphology, while the non co- 
planar merger, shown in the following Figure 16 and 17, 
is more extended and non disk-like. However, it shows 
the same time evolution of the surface brightness - the 
brightness in the IRAC bands increases as the brightness 
in the SCUBA bands decreases as more of the emitted 
energy shifts to shorter wavelengths. Finally, the set of 
images shown in Figures 18 and 19 depicts the time evo- 
lution of a simulation with a progenitor at z = 3, z3h270, 
(with a virial velocity comparable to h320 and e320, the 
ones shown in the previous two figures). While the trends 
in the time evolution of the surface brightness mentioned 
earlier for the z — scaled simulations are also seen here, 
one clear difference is that the z = 3 scaled simulations 
arc much more compact than the z = scaled simu- 
lations, which have emission extending to > 10 kpc at 
times, while the z = 3 simulation does not show any 
structure on scales larger than 1 kpc in the observed 
3.6 /zm band, though some extended long wavelength 
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emission is seen over several kpc in the observed 450 /im 
band. 

Chapman ct al. (2003a) obtained HST images of about 
a dozen SMGs (z ~ 2), and found that the morphologies 
were irregular and complex, suggestive of a merger ori- 
gin. Chapman et al. (2004) compared optical and radio 
imaging of this sample, and using the radio as a tracer of 
the star forming regions, argued for extended (about 10 
kpc, but not much larger), obscured starbursts in about 
70 % of their sample, while in the rest of the sample, 
the radio extends to about 1 kpc. While it is tempting 
to speculate about the differences in morphology between 
our z = and z = 3 mergers and the potential to witness 
such differences in morphology in high resolution images, 
it is not clear that we can use images alone (even of HST 
resolution) to determine whether our models for SMGs 
with disk properties of z = systems would be favored 
over those with properties of z = 3 systems. As we have 
noted earlier, a merger with progenitors representative of 
z = 3 galaxies in the Mo, Mao & White (1998) formal- 
ism, is less massive, even if it has the same virial velocity 
as a z = scaled merger. A more massive progenitor, 
even with disk properties of z — 3 simulations, can have a 
larger extent - hence, this mass degeneracy prohibits di- 
rect comparison between morphological structure seen in 
high resolution images and the the prevalence of mergers 
with progenitors at a given redshift. We note also that 
Genzel et al. (2003) studied a bright, lensed SMG at 
z = 2.8 to find a rotation velocity in excess of 400 km/s, 
and millimeter emission in a disk-like structure extend- 
ing to about 10 kpc. These properties seem to be akin 
to our h320, h500 simulations of co-planar mergers. It is 
important to note here that a disk-like image is a natural 
outcome of co-planar mergers (though such a orbital in- 
clination is expected to be rare); i.e., the appearance of 
a disk-like structure does not imply that the system had 
not undergone a major merger in the past. While there 
have been a number of recent papers (Pope et al. 2006; 
Iono et al. 2006) reporting high resolution imaging of 
SMGs, detailed modeling of individual sources requires 
kinematic information along the lines obtained by Genzel 
et al. (2006) for az~3 optically bright galaxy inferred 
to have a star formation rate of ~ 100 M Q /yr and a 
circular velocity of ~ 230 km/s. 

5. DISCUSSION 

In comparing our predictions to observed data, there 
are several important points of note. First, we have de- 
picted the IRAC and infrared X-ray correlations in the 
rest-frame, which leads to a clear dynamical interpreta- 
tion, as we have discussed previously. To test these pre- 
dictions will require a large number of observations in a 
narrow redshift slice. Second, in comparing the 850 /im 
fluxes (which have been shown at z = 2, unless otherwise 
noted) to observations, it will be necessary to disentangle 
the effects of lensing, and of multiple counterparts con- 
tributing to the wide SCUBA field 850 /im flux. More- 
over, a detailed study of current observational biases is 
needed. We defer this study to a future paper. Here, 
we point out a few basic observational comparisons - the 
IRAC color-color plot for galaxies placed at z = 0.3 — 2 
in the observed frame; the time-evolution of a rest-frame 
SED; a comparison of our SEDs to recent data obtained 
by Kovacs et al. (2006); the 850 /im flux (shown at z = 2) 



as function of the ratio of black hole luminosity to total 
luminosity, and as a function of the X-ray luminosity. 

Figure 20 shows the time evolution of the rest-frame 
SED for the A5e simulation. As the simulation pro- 
gresses, more of the energy is distributed to shorter wave- 
lengths, in a fashion similar to the time evolution of the 
SEDs for local ULIRGs as discussed by Chakrabarti et 
al. (2006a). Figures 14-19 depict this same trend in 
the surface brightness maps, as the apparent brightness 
in the 3.6 /im band increases as the brightness in the 
longer wavelength 450 /im band decreases. In Figure 21, 
we compare the SEDs from the simulations hl60, e226, 
e320, and e500 during the bright sub-mm phase, to ob- 
served SMG data obtained by Kovacs et al. (2006), and 
find reasonable agreement. In our simulations, except 
for mergers of very massive systems (Vv; r > 400 km/s), 
the brighter SMGs (F 850 ^ m > 1 mJy) at z ~ 2, do 
have non-negligible contribution from the black hole to 
the total bolometric luminosity. Kovacs et al. (2006) use 
the radio far-infrared correlation to find that the SMGs 
in their sample would not have a significant contribution 
from the black hole luminosity. Their sample is biased 
towards radio-loud sources, and does not include radio- 
quiet AGN that would contribute to the far-infrared lu- 
minosity, increasing the derived parameter, which is a 
measure of the far- infrared to radio luminosity. To better 
quantify our prediction for the variation of the 850 /im 
fluxes as a function of the black luminosity, we show in 
Figures 22 and 23 the 850 /im flux as a function of the 
ratio of the black hole luminosity to the total bolomet- 
ric luminosity and as a function of the X-ray luminosity. 
Figure 22 shows that there is a general trend for the 
850 /im flux to increase as a function of the ratio of the 
black hole luminosity to the total bolometric luminosity. 
The brighter SMGs which have low Z/BH/^totai are the 
more massive systems. A significant caveat to interpret- 
ing observations however, is a Compton thick population, 
which the current simulations cannot model in detail. 

We show in Figure 24 the IRAC color-color plot in 
the observed frame for simulations placed at z = 0.3, 1, 
and 2. The PAH spectral features redshift out of the 
IRAC bands by z <~ 0.5, and the relative increase in 
Ls /im/-^4.5 /im for energetically active AGN also can no 
longer be seen. However, the transition from a roughly 
spherical distribution of colors for low-redshift objects 
to a flattened elliptical distribution of colors for z <~ 2 
objects may be useful for interpreting observations. As 
mentioned previously, the z — 0.3 IRAC color-color plot 
displays the characteristic "bunny-ear" shape seen in the 
observed FLS color-color plot. We show the z = 0.3 
case here as Lacy et al. (2004) note that the median 
photometric redshift of their candidate AGN is ~ 0.3. 
This contrasts with the rest-frame IRAC color-color plot 
shown in Figure 11 primarily owing to the redshifting of 
the PAH template, which is a more sensitive function of 
redshift than the continuum. Ultimately, it will be use- 
ful to simulate observed color-color plots for a particular 
population which probe a range of redshifts. Work is un- 
derway to analyze current observational biases to realize 
this. 

Since the far-infrared (Lir, L70 /im) and to a lesser 
extent L 2 a /tin, are well correlated with the X-ray lumi- 
nosity, we review them here and give relations between 
these quantities. The red dots in Figure 25 identify cases 
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Fig. 14. — Photo Album of the Lifetime of an SMG produced in a co-planar merger (h320) in observed 3.6 (im band, (a) Pre-merger 
phase, (b) Close to Main Feedback Phase, (c) After Main Feedback Phase. 



Fig. 15. — Photo Album of the Lifetime of an SMG produced in a co-planar merger (h320) in observed 450 /im band at same times as in 
3.6 /im band. Times beyond the last phase would be black, i.e., would show no emission in SCUBA. This shows that more of the emitted 
energy is distributed to the shorter wavelengths as AGN feedback clears out the obscuring material. 



which have F^q > lmJy. The diagonal lines show 
constant Li/L x ; the best fit lines (the solid lines in these 
figures) are L 70 ^ m = 40 L x and L 2 4 (im = 10 L K respec- 
tively; there is nearly an order of magnitude of scatter 
relative to this line. 

The total integrated infrared luminosity and the hard 
X-ray luminosity are essentially global properties of these 
simulations. Hence, we investigate where SMGs from our 
merger simulations would lie on a Ltr — L K plane. Figure 
26 marks in red the location of SMGs (i.e., F 850 ^m > 
1 mjy for simulations placed at z = 2, at a luminosity 
distance of 15.5 Gpc) on the -Ltr — L x plane, with divi- 
sions into Class I, Class II, and Class III, which corre- 
spond to Ltr > 70 L x (the best-fit line that goes through 
our simulated data), Lm > 25 L x (the traditional quasar 
line, e.g., Elvis et al. 1994), and Lm < 5 i x respec- 
tively. Most of the SMGs studied by Alexander et al. 
(2005a) were found to have Lm ~ 200L X , on the ba- 
sis of the radio far-IR correlation. Indeed, quite a few 
from our simulations have similar _Ltr/Z/ x values, though 
most have lower values. Interestingly, samples of quasars 
(e.g. Page et al. 2001) have similar F 850 ^m values as 
SCUBA selected galaxies, which suggests that SMGs as 
defined by their 850 /im fluxes are a diverse and broad 
class of objects, as already highlighted by observational 



studies on the basis of their radio and optical properties 
(Ivison et al. 2000). Figure 26 demonstrates the diver- 
sity of SMGs, which traverse the Class I-Class II divide 
- some SMGs would fall on the traditional quasar line 
(£ir > 25I/ X ), yet many have higher values of Lir/L x . 
This shows explicitly that SMGs are a broader class of 
systems than quasars or starbursts (which have higher 
Lm/L K ), as selecting on the basis of the Fgso ^m values 
allows these objects to span a broad range in Ltr/ L x . We 
also show in Figure 26 a few observed sources from the 
literature (shown in yellow asterisks), comprised of local 
LIRGs (Mrk 463), ULIRGs (e.g., Mrk 231), quasars (e.g., 
3C273), and well-known SMGs (e.g., SMM J02399-0136) 
(as reported by Sanders & Mirabel 1996; Hughes et al. 
1997; Bautz et al. 2000; Genzel et al. 2003; Alexander 
et al. 2005a). Since we do not do a detailed comparison 
to observations in this paper, this plot certainly should 
not be interpreted to imply an exact correspondence be- 
tween our simulated data points and observed data in 
cases of overlap. Nonetheless, it is clear that z ~ 2 sys- 
tems such as SMM J02399-1036 (as studied in by Genzel 
et al. 2003) do have similar Ltr/L x ratios as the systems 
in our simulations, which suggests that the formation of 
such objects at z ~ 2 may be described by major merg- 
ers, whose time evolution is influenced by feedback from 
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FlG. 16. — Photo Album of the Lifetime of an SMG produced in a non co-planar merger (e320) in observed 3.6 ^.m band, (a) Pre-merger 
phase, (b) Close to Main Feedback Phase, (c) After Main Feedback Phase. 





Fig. 17. — Photo Album of the Lifetime of an SMG produced in a non co-planar merger (e320) in observed 450 fim band at same times 
as in shown above in the 3.6 fim band. Times beyond the last phase would be black; i.e., would show no emission in this band. 



the central AGN. As mentioned previously, the current 
merger simulations do not resolve the inner few parsecs, 
which may well be Compton thick and attenuate even the 
hard X-rays. However, Figure 26 does show that quasars 
like 3C273 in elliptical (McLeod & Rieke 1994b) host 
galaxies (though see recent work by Martel et al. 2003 
which identifies intricate features in the host galaxy not 
captured in early HST images), which are presumably 
more evolved, do lie closer to the Class III region than 
Mrk 463 or UGC 5101, which are interpreted to be in an 
earlier stage of a merger (Sanders & Mirabel 1996). 

A natural question to ask at this point is whether this 
classification scheme may be interpreted also as an evolu- 
tionary scheme within the context of the merger simula- 
tions studied here. Figure 27 presents the time evolution 
of the star formation rate, with the three classes marked, 
along with the time evolution of the 3.6 /jm image, and 
shows the track that a particular simulation (h320) would 
follow on the Lir — L x diagram (marked in yellow in the 
bottom panel of the Ltr/L x diagram). As shown, the 
system starts out at high values of Lir/L x ~ 100, and 
then after reaching a peak in L x , progressively moves 
down in Lir/L x towards the bottom portion of this dia- 
gram, which we have marked as Class III, and which we 
associate with elliptical merger remnants. (Other simula- 



tions evolve in a generally similar manner.) The Lir/L x 
values for ellipticals are not well determined, unfortu- 
nately. Extrapolating from SCUBA observations of a 
sample of elliptical galaxies (Di Matteo et al. 1999), we 
suggest that Lir ~ 5L X may represent an upper bound 
for ellipticals. We emphasize that the Class III designa- 
tion is highly uncertain owing to the lack of a determina- 
tion for Lir for merger remnants (though see Temi et al. 
2005 for a discussion of a correlation between the mid-IR 
emission and age of ellipticals) , as well as the wide range 
in X-ray luminosities (Mathews et al. 2006; O'Sullivan 
et al. 2001, 2004; Allen et al. 2000; Fabbiano 1989) and 
complex range of phenomena that produce X-ray lumi- 
nosities in ellipticals, i.e., X-ray emission from hot gas, 
black hole, or residual star formation. Nonetheless, the 
increase in L x /Lb with age (Mackie & Fabbiano 1997; 
O'Sullivan et al. 2001) does suggest that the amount 
of X-ray emitting hot gas increases with time relative 
to the amount of cold (infrared emitting) gas. It will 
also be useful to directly compare correlations between 
photometric and kinematic correlations for merger rem- 
nants, as explored recently by Rothberg & Joseph (2006) 
to better quantify this tentative evolutionary scheme. 
We study correlations between kinematic and multiwave- 
length photometric properties in a future paper, and an- 
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Fig. 18. — Photo Album of the Lifetime of an SMG produced of a z = 3 scaled merger (z3e226) in observed 3.6 /an band, (a) Pre-merger 
phase, (b) Close to Main Feedback Phase, (c) After Main Feedback Phase. 



Fig. 19. — Photo Album of the Lifetime of an SMG produced of a z = 3 scaled merger (z3e226) in observed 450 (im band. Times beyond 
the last phase would be black; i.e., would show no emission in this band. 




Fig. 20.— Time Evolution of a Rest-Frame SED for e320. As 
the AGN feedback clears out the obscuring material, more of the 
emitted energy shifts from the longer wavelengths to the shorter 
wavelengths (compare the solid line which is the pre-merger phase 
to the dash-dotted line which is from the post-feedback phase). 



Fig. 21. — Comparison to SHARC-2 and SCUBA data reported 
in Kovacs et al., astro-ph/0604591 Dash-double-dotted line is the 
hl60 simulation, solid line is e22t>, dashed line is e320, and dash- 
dotted line is e500; all are shown close to the peak of their sub-mm 
bright phase. 
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Fig. 22. — The observed Fgso fim flux ( m Jy) as a function of the 
ratio of black hole luminosity to total luminosity. 
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Fig. 23. — The observed Fg^o p, m flux (mJy) as a function of the 
X-ray luminosity. 



alyze a larger set of simulations to quantitatively address 
the transition from SMGs to ellipticals. From our pre- 
liminary analysis here, Figure 27 does suggest that the 
evolution of SMGs, from the pre-merger phase to the 
merger remnant, may be qualitatively understood in a 
relatively simple manner, wherein the system occupies 
some characteristic region of the L\n — L K diagram dur- 
ing each phase of its lifetime. 

6. CONCLUSION 

• Our simulations of gas-rich major mergers with 
black hole feedback naturally lead to the production 
of SMGs, which match photometric observations. The 
SMGs formed in these simulations have star formation 
rates of ~ 500 — 1000 M Q /yr, infrared luminosities of 
~ l-5xlO 12 i , and virial velocities of ~ 300— 400km/s. 

• We comment on the M sta r — Mbh relation for SMGs, 



and its inference from observations. We do find that the 
black holes in SMGs are in a rapidly growing phase, and 
grow by factors ~ 5 between z ~ 2 and the present day. 
However, we do not find that SMGs at redshift 2 would 
lie two orders of magnitude below the local relation which 
is an extreme lower limit that follows from deriving the 
black hole masses under the assumption of Eddington- 
limited accretion. 

• We demonstrate that clustering in IRAC color-color 
space can be naturally explained within the context of 
the merger simulations studied here; clustering in this 
context translates to the system spending more of its 
lifetime in a given region of color-color space; it is also 
positively correlated with the stars dominating in their 
contribution to the bolometric luminosity. We recover 
a similar percentage of sources occupying the AGN- 
demarcated region as in the Lacy et al. (2004) plot. We 
use a phenomenological model for the PAH emission to 
demonstrate the change in color-color plot when a PAH 
template is included in calculating the colors, both in the 
rest-frame, and in the observed frame for galaxies from 
z = 0.3-2. 

• Our models would predict an inherent similarity 
between color (not fluxes) evolution for local ULIRGs 
and high redshift luminous galaxies, which is established 
largely by the dynamical effect of feedback from the AGN 
in these simulations. The extension of this model to high 
rcdshifts (z > 6) would suggest that these correlations 
would also hold for the highest redshift galaxies known. 

• We predict a correlation between the rest-frame in- 
frared, 70 fim, and 24 /im luminosity and the hard X- 
ray luminosity for SMGs, with an increase in scatter at 
low X-ray luminosities. To aid observational studies, we 
quantify the far-infrared X-ray correlations. These pre- 
dictions will be directly testable by future instruments, 
such as Herschel, and possibly through source stacking 
analysis using current Spitzer data. 

• Our photo albums of the lifetimes of SMGs visually 
illustrate the differential variation in surface brightness 
between the SCUBA and IRAC bands. Of particular 
note is the increase in apparent brightness in the IRAC 
bands, which is concomitant with a decrease in bright- 
ness in the SCUBA bands, towards the late phases of the 
merger. This shows that more and more of the emitted 
energy shifts to the shorter wavelengths as AGN feedback 
disperses the obscuring material. We also demonstrate 
that the morphology as seen in these bands is partly a 
function of orbital inclination, with co-planar mergers 
producing disk-like morphologies in the active phase. As 
such, even galaxies with observed disk-like morphologies 
may have experienced major (co-planar) mergers. The 
simulations of merger with progenitors at high rcdshifts 
lead to more compact morphologies than progenitors of 
lower redshift systems, which generally have emission ex- 
tending to ~ 10 kpc in the active phase in the IRAC and 
SCUBA bands. 

• We find that our predicted SEDs are in good agree- 
ment with recent multiwavelength photometry. We show 
the IRAC color-color plot for galaxies in a narrow red- 
shift slice, with z ~ 2, which results in a flattened el- 
liptical distribution of colors unlike the nearly spherical 
distribution at low redshifts. 

• We depict the variation of the observed 850 /xm flux 
as a function of both the X-ray luminosity and the ratio 
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Fig. 24. — IRAC color-color plot in observed frame, with all simulations placed at a) z = 0.3, b) z = 1, c) z = 2. Crosses designate the 
continuum and squares the combined spectrum, i.e., the continuum and the PAH model, as before. Note the slightly "bunny-ear" shape of 
the IRAC color-color plot for the z = 0.3 slice. 




Fig. 25. — (a) L70 (rest-frame) vs. the hard X-ray luminosity, and (b) L24 fim (rest-frame) vs. the hard X-ray luminosity. Red dots 
denote cases which have i"850 /an > lmjy. The solid lines correspond to L70 fim = 40 L x and L24 (im = 10 L x . 



of the black hole luminosity to the total bolometric lumi- 
nosity. There is a general trend for SMGs on the bright 
end at z ~ 2 to have a significant (> 50%) contribution 
from the black hole to their total bolometric luminosity. 

• We stress that the correlations and photometric prop- 
erties of SMGs listed here are predictions of our model 
and can be observationally tested - by obtaining a large 
sample (~ 100) of observations in a narrow rcdshift 
range. 

• We find that SMGs are a broader class of systems 
than quasars or starbursts. We introduce a simple, 
heuristic classification scheme for SMGs on the basis of 
their Lir/L x ratios, Lm > 70 L x is Class I, Ltr > 25 L x 
is Class II, Liji < 5 L x is Class III. We suggest that 
this may also be interpreted as an evolutionary scheme 
as SMGs transit from the pre-merger stage through the 



quasar phase to merger remnants. 
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Fig. 27. — An Evolutionary Model for SMGs: (a) top panel shows the time evolution of the star formation rate, with the Class I — > Class 
II — > Class III time evolution marked. Red dots mark the times of the images shown below, (b) middle panel shows the 3.6 jim images 
of h320 along this time sequence, (c) bottom panel shows the time evolution on the Lir — L x plot in yellow for this simulation. As time 
progress, Ltr/L x increases first as the galaxies merge and the black hole becomes energetically active, and decreases as black hole feedback 
expels the obscuring material, ultimately producing a merger remnant. 
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